
* Do-file to replicate tables in "Estimating the demand for business training: Evidence from Jamaica"


cap log close
clear all
set more off

* Set directory
cap cd {INSERT PATH HERE}

* Set global for subdirectories
//global main_dir c(pwd)
global data_raw "Data/Raw data"
global data_int "Data/Intermediate data"
global data_fin "Data/Final data"
global results_tab "Results/Tables"
global results_fig "Results/Figures"
global codes "Codes"


* Set global options for graphs (font if Stata >15 and style)
cap graph set window fontface "Times New Roman"



*****************************************************
*		FIGURES
*****************************************************

***********************************************************
*** Figure 1: Demand curve, WTP vs payments ***
***********************************************************

use "$data_fin/data_anonymized", clear
drop if price==0
replace price=price/1000
replace wtp_real=wtp_real/1000

gen pricet=price if tioli_spl==1
gen priceb=price if tioli_spl==0


* Generate mean of payment variables (instead of collapsing the dataset)
bys pricet: egen mean_paid_fullt=mean(paid_full) if tioli_spl==1
bys pricet: replace mean_paid_fullt=. if _n!=1
bys priceb: egen mean_paid_fullb=mean(paid_full) if tioli_spl==0
bys priceb: replace mean_paid_fullb=. if _n!=1

* Generate share of people willing to buy
cumul wtp_real, gen(cumwtp)
gen share=1-cumwtp
sort share		

* Plot demand curve comparing stated WTP and actual payment
twoway(line mean_paid_fullb priceb if inrange(price,1,20) & tioli_spl==0, sort lpattern(longdash)) (line mean_paid_fullt pricet if tioli_spl==1 & inrange(price,5,15), sort lpattern(shortdash)) (line share wtp_real if inrange(wtp_real,1.5,30), sort) , ///
	ylab(, grid) ytitle("Proportion of firms", axis(1))  ///
	xlabel(1 5(5)30, format(%15.0fc)) xtitle("Price in J$000") ///
	xline(20, lpattern(solid)) text(0.4 20 "Cost of Course", placement(ne)) ///
	legend(order(1 "Paid (BDM)" 2 "Paid (TIOLI)" 3 "WTP (BDM)")) graphregion(fcolor(white))
graph export "$results_fig/figure1n.pdf", replace



***********************************************************
*** Figure 2: WTP and Demand by Treatment               ***
***********************************************************

*****************************************
* Panel A. WTP by Treatment: BDM Sample
*****************************************
use "$data_fin/data_anonymized", clear
keep if tioli_spl==0
replace wtp_real=wtp_real/1000

*winsorize only 2 cases for the graph
replace wtp_real=25 if wtp_real>25

* Generate share of people willing to buy by treatment group
cumul wtp_real if form==1, gen(cumwtp1) equal
gen share1=1-cumwtp1
cumul wtp_real if form==2, gen(cumwtp2) equal
gen share2=1-cumwtp2
cumul wtp_real if form==3, gen(cumwtp3) equal
gen share3=1-cumwtp3
cumul wtp_real if form==4, gen(cumwtp4) equal
gen share4=1-cumwtp4
sort wtp_real

* Plot demand curve by treatment arm
twoway (line share1 wtp_real if form==1) (line share2 wtp_real if form==2, lpattern(-)) (line share3 wtp_real if form==3, lpattern(_)) (line share4 wtp_real if form==4, lpattern(-.-)), ///
	ylab(, grid) ytitle("Share willing to purchase") title("WTP by treatment: BDM sample") ///
	xlabel(0(5)25, format(%15.0fc)) xtitle("Willingness to Pay in J$000") ///
	xline(20, lpattern(dash)) text(0.4 20 "Cost of Course", placement(ne)) ///
	legend(order(1 "Basic" 2 "Credit" 3 "Discount Frame" 4 "Credit and Discount Frame")) graphregion(fcolor(white)) 
graph save "$results_fig/bdm_demand_bytreat.gph", replace


*****************************************
* Panel B. Demand by Treatment
*****************************************

use "$data_fin/data_anonymized", clear
drop if price==0
replace price=price/1000

* Generate unique credit dummy (for both BDM and TIOLI)
gen credit_unique=(form==2 | form==4 | credit==1)

* Collapse data 
collapse paid_full paid_something, by(tioli_spl price credit_unique)


* Plot actual payments by credit (superimposed)
twoway (line paid_full price if tioli_spl==0 & credit==0, sort lpattern(dash)) (line paid_full price if tioli_spl==0 & credit==1, sort) ///
	(line paid_full price if tioli_spl==1 & credit==0, sort lpattern(dash)) (line paid_full price if tioli_spl==1 & credit==1, sort), ///
	ylab(, grid) ytitle("Share that paid full price") ///
	xlabel(1 3 5 10 15 20, format(%15.0fc)) xtitle("Price in J$") ///
	legend(label(1 "No credit, BDM sample") label(2 "Credit, BDM sample") label(3 "No credit, TIOLI sample") label(4 "Credit, TIOLI sample")) ///
	graphregion(fcolor(white))
graph save "$results_fig/all_demand_bycredit_overlap.gph", replace

* Plot actual payments by credit (side-by-side)
twoway (line paid_full price if tioli_spl==0 & credit==0, sort lpattern(dash)) (line paid_full price if tioli_spl==0 & credit==1, sort), ///
	ylab(, grid) ytitle("Share that paid full price") ///
	xlabel(1 3 5 10 15 20, format(%15.0fc)) xtitle("Price in J$000") ///
	title("Demand: BDM sample") legend(label(1 "No credit") label(2 "Credit") rows(2)) ///
	graphregion(fcolor(white)) name(bdm, replace)
twoway (line paid_full price if tioli_spl==1 & credit==0, sort lpattern(dash)) (line paid_full price if tioli_spl==1 & credit==1, sort), ///
	ylab(, grid) ytitle("Share that paid full price") ///
	xlabel(5 10 15, format(%15.0fc)) xtitle("Price in J$000") ///
	title("Demand: TIOLI sample") legend(label(1 "No credit") label(2 "Credit") rows(2)) ///
	graphregion(fcolor(white)) name(tio, replace)
graph combine bdm tio, graphregion(margin(zero)) xsize(8)
graph save "$results_fig/all_demand_bycredit_byside.gph", replace


*Combine graphs
graph combine "$results_fig/bdm_demand_bytreat.gph" "$results_fig/all_demand_bycredit_byside.gph"  , rows(2)   graphregion(fcolor(white))
graph export "$results_fig/FigureA1n.pdf", replace


********************************************************
*** Figure 3: Demand by Covariates  BDM               **
********************************************************

use "$data_fin/data_anonymized", clear
keep if tioli_spl==0
*Gender
replace price=price/1000
drop if price==0
collapse paid_full, by(price female)
twoway (line paid_full price if female==1, lpattern(-)) (line paid_full price if female==0 ), ///
	ylab(, grid) ytitle("Share purchasing") ///
	xlabel(1 3 5(5)20, format(%15.0fc)) xtitle("Price in J$000") ///
	  text(.35 15 "p = 0.989", place(e) size(.4cm)) ///
	legend(order(1 "Female" 2 "Male") rows(2)) graphregion(fcolor(white)) saving("$results_fig/Figure2anp.gph", replace) title("Demand by Gender")

	
*Wealth
use "$data_fin/data_anonymized", clear
keep if tioli_spl==0
replace price=price/1000
drop if price==0
sum wealthindex,det
gen pwealth=(wealthindex>=r(p50)) if wealthindex!=.
collapse paid_full, by( price pwealth )
twoway (line paid_full price if pwealth==0, lpattern(-)) (line paid_full price if pwealth==1), ///
	ylab(, grid) ytitle("Share purchasing") ///
	xlabel(1 3 5(5)20, format(%15.0fc)) xtitle("Price in J$000") ///
	text(.35 15 "p = 0.446", place(e) size(.4cm)) ///
	legend(order(1 "< median wealth" 2 "> median wealth") rows(2)) graphregion(fcolor(white)) saving("$results_fig/Figure2bnp.gph", replace) title("Demand by Wealth")


* Sales
use "$data_fin/data_anonymized", clear
keep if tioli_spl==0
replace price=price/1000
drop if price==0
drop if sales_dollars==.
sum sales_dollars, det
gen lowsales=sales_dollars<r(p50)
gen highsales=sales_dollars>=r(p50)
collapse paid_full, by( price lowsales highsales )

	twoway (line  paid_full price if lowsales==1 , lpattern(-)) (line  paid_full price if highsales==1 ), ///
 ylab(, grid) ytitle("Share purchasing") ///
	xlabel(1 3 5(5)20, format(%15.0fc)) xtitle("Price in J$000") ///
	 text(.35 15 "p = 0.042", place(e) size(.4cm)) ///
		legend(order(1  "< median sales" 2 "> median sales") rows(3)) graphregion(fcolor(white)) saving("$results_fig/Figure2dnp.gph", replace) title("Demand by Firm Size")

* Education
use "$data_fin/data_anonymized", clear
keep if tioli_spl==0
replace price=price/1000
drop if price==0
collapse paid_full, by( price edu_comp )

twoway (line paid_full price  if edu_comp==0 , lpattern(-))  (line paid_full price  if edu_comp==1  ), ///
	ylab(, grid) ytitle("Share purchasing") ///
	xlabel(1 3 5(5)20, format(%15.0fc)) xtitle("Price in J$000") ///
	 text(.35 15 "p = 0.001", place(e) size(.4cm)) ///
	legend(order(1 "Low Education" 2 "High Education") rows(2)) graphregion(fcolor(white)) saving("$results_fig/Figure2gnp.gph", replace) title("Demand by Education")
	
graph combine "$results_fig/Figure2anp.gph" "$results_fig/Figure2bnp.gph"  "$results_fig/Figure2dnp.gph"  "$results_fig/Figure2gnp.gph" , col(2) graphregion(fcolor(white))
graph export "$results_fig/Figure2np.pdf", replace


********************************************************
*** Figure 4: Demand by Covariates  TIOLI             **
********************************************************

use "$data_fin/data_anonymized", clear
keep if tioli_spl==1
*Gender
replace price=price/1000
drop if price==0
collapse paid_full, by(price female)
twoway (line paid_full price if female==1, lpattern(-)) (line paid_full price if female==0 ), ///
	ylab(, grid) ytitle("Share purchasing") ///
	xlabel(5(5)15, format(%15.0fc)) xtitle("Price in J$000") ///
	  text(.50 10 "p = 0.167", place(e) size(.4cm)) ///
	legend(order(1 "Female" 2 "Male") rows(2)) graphregion(fcolor(white)) saving("$results_fig/Figure2anpt.gph", replace) title("Demand by Gender")

	
*Wealth
use "$data_fin/data_anonymized", clear
keep if tioli_spl==1
replace price=price/1000
drop if price==0
sum wealthindex,det
gen pwealth=(wealthindex>=r(p50)) if wealthindex!=.
collapse paid_full, by( price pwealth )
twoway (line paid_full price if pwealth==0, lpattern(-)) (line paid_full price if pwealth==1), ///
	ylab(, grid) ytitle("Share purchasing") ///
	xlabel(5(5)15, format(%15.0fc)) xtitle("Price in J$000") ///
	text(.50 10 "p = 0.004", place(e) size(.4cm)) ///
	legend(order(1 "< median wealth" 2 "> median wealth") rows(2)) graphregion(fcolor(white)) saving("$results_fig/Figure2bnpt.gph", replace) title("Demand by Wealth")


* Sales
use "$data_fin/data_anonymized", clear
keep if tioli_spl==1
replace price=price/1000
drop if price==0
drop if sales_dollars==.
sum sales_dollars, det
gen lowsales=sales_dollars<r(p50)
gen highsales=sales_dollars>=r(p50)
collapse paid_full, by( price lowsales highsales )

	twoway (line  paid_full price if lowsales==1 , lpattern(-)) (line  paid_full price if highsales==1 ), ///
 ylab(, grid) ytitle("Share purchasing") ///
	xlabel(5(5)15, format(%15.0fc)) xtitle("Price in J$000") ///
	 text(.50 10 "p = 0.586", place(e) size(.4cm)) ///
		legend(order(1  "< median sales" 2 "> median sales") rows(3)) graphregion(fcolor(white)) saving("$results_fig/Figure2dnpt.gph", replace) title("Demand by Firm Size")

* Education
use "$data_fin/data_anonymized", clear
keep if tioli_spl==1
replace price=price/1000
drop if price==0
collapse paid_full, by(price edu_comp )

twoway (line paid_full price  if edu_comp==0 , lpattern(-))  (line paid_full price  if edu_comp==1  ), ///
	ylab(, grid) ytitle("Share purchasing") ///
	xlabel(5(5)15, format(%15.0fc)) xtitle("Price in J$000") ///
	 text(.50 10 "p = 0.157", place(e) size(.4cm)) ///
	legend(order(1 "Low Education" 2 "High Education") rows(2)) graphregion(fcolor(white)) saving("$results_fig/Figure2gnpt.gph", replace) title("Demand by Education")
	
graph combine "$results_fig/Figure2anpt.gph" "$results_fig/Figure2bnpt.gph"  "$results_fig/Figure2dnpt.gph"  "$results_fig/Figure2gnpt.gph" , col(2) graphregion(fcolor(white))
graph export "$results_fig/Figure2npt.pdf", replace


*********************************************		
**** Figure 5 Sunk Costs. BDM Sample      ***
*********************************************

use "$data_fin/data_anonymized", clear
keep if tioli_spl==0

* Rescale price 
gen price_k=price/1000
lab var price_k "Price (000)"

* Relabel original price var
lab def price 0 "Price = 0" 1000 "Price = 1,000" 3000 "Price = 3,000" 5000 "Price = 5,000" 10000 "Price = 10,000" 15000 "Price = 15,000" 20000 "Price = 20,000", modify

*Number of classes. Pooled Sample
preserve
keep if paid_full==1
collapse (mean) attendance_total  , by(price_k)
sort price_k
twoway line attendance_total price_k, xtitle("Price in J$000") ytitle("Number of classes attended") ///
	legend(order(1 "No controls")) graphregion(fcolor(white)) 
restore
graph save "$results_fig/fig3na.gph", replace

preserve
keep if paid_full==1
collapse (mean) atleast1class , by(price_k)
sort price_k
twoway line atleast1class price_k, xtitle("Price in J$000") ytitle("At least 1 class attended") ///
	legend(order(1 "No controls")) graphregion(fcolor(white)) 
restore
graph save "$results_fig/fig3nb.gph", replace

graph combine  "$results_fig/fig3nb.gph" "$results_fig/fig3na.gph"  , col(2) graphregion(fcolor(white))
graph export "$results_fig/figure3n.pdf", replace



**********************************
* Figures A1-A3
**********************************

use "$data_fin/data_anonymized", clear
replace price=price/1000
replace wtp_real=wtp_real/1000

gen bdm_spl=1-tioli_spl


*Dummy for any missing covariate that will be used in matching
gen mis_cov=0
foreach var of varlist female wealthindex sales_lm  edu_compsecondary risk_taking age  expect_inc1020 expect_inc20 pi_index practice_share   {
	replace mis_cov=1 if `var'==.
}
*Replacing missing with mean value or 1 for dummies, only for matching
cap drop m_* d_*
* Compute variables to replace missing values 
foreach var in  wealthindex sales_lm  risk_taking age  pi_index practice_share {
		clonevar mat_`var'=`var'
		sum  mat_`var'
		replace mat_`var'=r(mean) if missing(`var')
	}
foreach var in  expect_inc1020 expect_inc20  edu_compsecondary {
		clonevar mat_`var'=`var'
		replace mat_`var'=1 if missing(`var')
	}


label var mat_age "Age"
label var mat_edu_compsecondary "Education: more than secondary"
label var mat_wealthindex "Wealth Index"
label var mat_risk_taking "Risk taking index [0-10]"
label var  mat_pi_index "Personal initiative index [1-7]"
label var mat_sales_lm "Sales last month"
label var mat_practice_share "Share of bus. practices"
label var mat_expect_inc1020  "Exp. sales incr. 10%-20%"
label var mat_expect_inc20 "Exp. sales incr. by more than 20%"
label var mis_cov "Missing values in covariates"


*Estimating ATT using propensity score kernel matching, logit model for pscore
*Treatment variable: BDM sample, outcome: paid full, we match to make BDM sample and TIOLI sample similar
drop if price==0
gen pricet=price if tioli_spl==1
gen priceb=price if tioli_spl==0

*Using only wealthindex
gen mis_wealth=wealthindex==.
kmatch ps bdm_spl  mat_wealthindex mis_wealth (paid_full) if price>=5 & price<=15, att po vce(boot) over(price) comsup
mat Ew=e(b)

gen mean_paid_wealthbdm=.
gen mean_paid_wealthtioli=.
replace mean_paid_wealthbdm=matrix(Ew[1,2]) if price==5
replace mean_paid_wealthtioli=matrix(Ew[1,3]) if price==5
replace mean_paid_wealthbdm=matrix(Ew[1,5]) if price==10
replace mean_paid_wealthtioli=matrix(Ew[1,6]) if price==10
replace mean_paid_wealthbdm=matrix(Ew[1,8]) if price==15
replace mean_paid_wealthtioli=matrix(Ew[1,9]) if price==15
bys priceb: replace mean_paid_wealthbdm=. if _n!=1
bys pricet: replace mean_paid_wealthtioli=. if _n!=1


*Using only age
gen mis_age=age==.
kmatch ps bdm_spl mat_age mis_age (paid_full) if price>=5 & price<=15, att po vce(boot) over(price) comsup
mat Ea=e(b)
gen mean_paid_agebdm=.
gen mean_paid_agetioli=.
replace mean_paid_agebdm=matrix(Ea[1,2]) if price==5
replace mean_paid_agetioli=matrix(Ea[1,3]) if price==5
replace mean_paid_agebdm=matrix(Ea[1,5]) if price==10
replace mean_paid_agetioli=matrix(Ea[1,6]) if price==10
replace mean_paid_agebdm=matrix(Ea[1,8]) if price==15
replace mean_paid_agetioli=matrix(Ea[1,9]) if price==15
bys priceb: replace mean_paid_agebdm=. if _n!=1
bys pricet: replace mean_paid_agetioli=. if _n!=1

*Using only education, and exact matching for the dummy on complete highschool
gen mis_edu=edu_compsecondary==.
kmatch em bdm_spl mat_edu_compsecondary mis_edu (paid_full) if price>=5 & price<=15, att po vce(boot) over(price) 
mat Ee=e(b)
gen mean_paid_edubdm=.
gen mean_paid_edutioli=.
replace mean_paid_edubdm=matrix(Ee[1,2]) if price==5
replace mean_paid_edutioli=matrix(Ee[1,3]) if price==5
replace mean_paid_edubdm=matrix(Ee[1,5]) if price==10
replace mean_paid_edutioli=matrix(Ee[1,6]) if price==10
replace mean_paid_edubdm=matrix(Ee[1,8]) if price==15
replace mean_paid_edutioli=matrix(Ee[1,9]) if price==15
bys priceb: replace mean_paid_edubdm=. if _n!=1
bys pricet: replace mean_paid_edutioli=. if _n!=1

*Using all variables in the demand table for the match
kmatch ps bdm_spl female mat_wealthindex mat_sales_lm mat_risk_taking mat_age mat_pi_index mat_practice_share mat_expect_inc1020 mat_expect_inc20 mat_edu_compsecondary mis_cov (paid_full) if price>=5 & price<=15, att po vce(boot) over(price) comsup 
mat Ef=e(b)
gen mean_paid_fullbdm=.
gen mean_paid_fulltioli=.
replace mean_paid_fullbdm=matrix(Ef[1,2]) if price==5
replace mean_paid_fulltioli=matrix(Ef[1,3]) if price==5
replace mean_paid_fullbdm=matrix(Ef[1,5]) if price==10
replace mean_paid_fulltioli=matrix(Ef[1,6]) if price==10
replace mean_paid_fullbdm=matrix(Ef[1,8]) if price==15
replace mean_paid_fulltioli=matrix(Ef[1,9]) if price==15
bys priceb: replace mean_paid_fullbdm=. if _n!=1
bys pricet: replace mean_paid_fulltioli=. if _n!=1
kmatch summarize
mat M = r(M)

*Plot pooling differences in covariates across matched and non-matched for the 3 prices
coefplot (matrix(M[,3]), msymbol(triangle)) (matrix(M[,6]) , msymbol(circle)) , legend(order(2 "Raw" 4 "Matched")) xline(0) xtitle(Std. mean difference)  graphregion(fcolor(white))
graph save "$results_fig/pscore_matching_stats.gph", replace
graph export "$results_fig/pscore_matching_stats.pdf", replace
kmatch density, lw(*6 *2) lc(*.5 *1)  graphregion(fcolor(white))
graph save  "$results_fig/pscore_matching_dist.gph", replace
graph export  "$results_fig/pscore_matching_dist.pdf", replace
*Note: we manually changed background to white color (not allowed by command) and saved as  "$results_fig/pscore_matching_distw.pdf"

* Generate mean of payment variables
bys pricet: egen mean_paid_fullt=mean(paid_full) if tioli_spl==1
bys pricet: replace mean_paid_fullt=. if _n!=1
bys priceb: egen mean_paid_fullb=mean(paid_full) if tioli_spl==0
bys priceb: replace mean_paid_fullb=. if _n!=1


foreach var in age wealth edu full  {

* Plot demand curves and re-weigthed curves using potential outcomes from matching algorithm
twoway(line mean_paid_fullb price if inrange(price,1,20) & tioli_spl==0, sort lpattern(longdash)) (line mean_paid_fullt pricet if tioli_spl==1 & inrange(price,5,15), sort lpattern(shortdash)) (line mean_paid_`var'bdm priceb if inrange(price,5,15) & tioli_spl==0, sort lpattern(solid)) (line mean_paid_`var'tioli pricet if inrange(price,5,15) & tioli_spl==1, sort lpattern(dash)) , ///
	ylab(, grid) ytitle("Proportion of firms", axis(1))  ///
	xlabel(1 5(5)30, format(%15.0fc)) xtitle("Price in J$000") ///
	xline(20, lpattern(solid)) text(0.4 20 "Cost of Course", placement(ne)) ///
	legend(order(1 "BDM" 2 "TIOLI" 3 "BDM `var'" 4 "TIOLI `var'" )) graphregion(fcolor(white))
graph save "$results_fig/figure1ex`var'.gph", replace
}

graph combine "$results_fig/Figure1exage.gph" "$results_fig/Figure1exedu.gph"  "$results_fig/Figure1exwealth.gph"  "$results_fig/Figure1exfull.gph",   col(2) graphregion(fcolor(white))
graph export "$results_fig/Figure1ex.pdf", replace





*********************************************************************
*** 					TABLES                                    ***
********************************************************************* 

*******************************************************
*** Table 1: Balance by price offered               ***
*******************************************************
use "$data_fin/data_anonymized", clear

* Set globals with lists of covariates for balance table
global spec1 expect_inc1020 expect_inc20 finance20k_notpos
global spec2 practice_share has_employees msect_services msect_retail msect_agric vehicle_bus internet
global spec3 age female married edu_compsecondary pi_index risk_taking training_previous
global spec4 practices_comparison proactive_comparison res_wage wtp

global nspec1  female age married edu_compsecondary  internet wealthindex  
global nspec2  training_previous wtp risk_taking  pi_index  practices_comparison proactive_comparison
global nspec3  has_employees msect_services msect_retail msect_manuf sales_lm profits_lm
global nspec4  practice_share expect_inc1020 expect_inc20 finance20k_notpos

* Set globals with lists of covariates for summary stats table
global vars1 complete_online active_business busage_0to1 busage_0 busage_1to5 busage_5more ///
	has_buspartner has_employees employees_total whours_worked sector* records_keeps practice_share sales_positive sales_dollars ///
	profits_negative profits_zero profits_positive profits_dollars owns_othercomp finance20k* vehicle_bus
global vars2 female age edu_compsecondary married internet res_wage_dollars subsist_entrep pi_index risk_taking ///
	practices_comparison proactive_comparison expect* wtp_dollars training_previous ///
	mngbus* loc1-loc7

global f_vars1 f_selfemp f_has_employees f_employees_total f_whours_worked f_practice_share f_sales_poslm f_sales_lastm_dollars ///
	f_profits_neglm f_profits_zerolm f_profits_poslm f_profits_lastm_dollars f_owns_othercomp
global f_vars2 f_pi_index
global f_vars3 f_know_1 f_know_2 f_know_3 f_know_4 f_know_5 f_knowledge
global f_vars4 f_course_took f_course_satisfac f_course_payment f_course_use f_course_metppl f_course_wtp f_course_fairprice


cap drop m_* d_*
* Compute variables to replace missing values (for joint test)
global m_varlist ""
global d_varlist ""
foreach var in $nspec1 $nspec2 $nspec3 $nspec4  {
	clonevar m_`var'=`var'
	replace m_`var'=0 if missing(`var')
	gen d_`var'=missing(`var')
	global m_varlist "$m_varlist m_`var'"
	global d_varlist "$d_varlist d_`var'"
}

* Rescale price not to have always 0 coefficients
replace price=price/1000

eststo clear

* Store summary statistics for BDM sample
eststo bdm: estpost summ $nspec1 $nspec2 $nspec3 $nspec4 if tioli_spl==0

* Save results of regressions in a matrix and add to stored results for BDM sample
tempname mat_diff mat_se mat_pval
foreach var in $nspec1 $nspec2 $nspec3 $nspec4  {
	qui reg `var' price if tioli_spl==0, robust
	matrix `mat_diff'=nullmat(`mat_diff'),_b[price]			// store beta (difference)
	matrix `mat_se'=nullmat(`mat_se'),_se[price]			// store standard error
	qui test price
	matrix `mat_pval'=nullmat(`mat_pval'),r(p)				// store p-value (for the sign. stars)
}
matrix colnames `mat_diff' = $nspec1 $nspec2 $nspec3 $nspec4 
matrix colnames `mat_se' = $nspec1 $nspec2 $nspec3 $nspec4 
matrix colnames `mat_pval' = $nspec1 $nspec2 $nspec3 $nspec4 
estadd matrix diff=`mat_diff': bdm
estadd matrix se=`mat_se': bdm
estadd matrix pval=`mat_pval': bdm
matrix drop `mat_diff' `mat_se' `mat_pval'

* Add number of observations for BDM sample
qui count if tioli_spl==0
estadd scalar N=r(N), replace: bdm


* Store summary statistics for TIOLI sample
eststo tioli: estpost summ $nspec1 $nspec2 $nspec3 $nspec4  if tioli_spl==1

* Save results of regressions in a matrix and add to stored results for TIOLI sample
tempname mat_diff mat_se mat_pval
foreach var in $nspec1 $nspec2 $nspec3 $nspec4  {
	qui reg `var' price i.demo_session if tioli_spl==1, robust
	matrix `mat_diff'=nullmat(`mat_diff'),_b[price]			// store beta (difference)
	matrix `mat_se'=nullmat(`mat_se'),_se[price]			// store standard errors
	qui test price
	matrix `mat_pval'=nullmat(`mat_pval'),r(p)				// store p-value (for the sign. stars)
}
matrix colnames `mat_diff' = $nspec1 $nspec2 $nspec3 $nspec4 
matrix colnames `mat_se' = $nspec1 $nspec2 $nspec3 $nspec4 
matrix colnames `mat_pval' = $nspec1 $nspec2 $nspec3 $nspec4 
estadd matrix diff=`mat_diff': tioli
estadd matrix se=`mat_se': tioli
estadd matrix pval=`mat_pval': tioli
matrix drop `mat_diff' `mat_se' `mat_pval'

* Add number of observations for TIOLI sample
qui count if tioli_spl==1
estadd scalar N=r(N), replace: tioli


* Test differences between BSM and TIOLI and add to new matrix
eststo bdm_tioli: estpost summ $nspec1 $nspec2 $nspec3 $nspec4 	// these results are not printed out, they are just to have stored results to append the matrices below
tempname mat_diff mat_se mat_pval
foreach var in $nspec1 $nspec2 $nspec3 $nspec4  {
	qui reg `var' tioli_spl, robust
	matrix `mat_diff'=nullmat(`mat_diff'),_b[tioli_spl]			// store beta (difference)
	matrix `mat_se'=nullmat(`mat_se'),_se[tioli_spl]			// store standard errors
	qui test tioli_spl
	matrix `mat_pval'=nullmat(`mat_pval'),r(p)					// store p-value (for the sign. stars)
}
matrix colnames `mat_diff' = $nspec1 $nspec2 $nspec3 $nspec4 
matrix colnames `mat_se' = $nspec1 $nspec2 $nspec3 $nspec4 
matrix colnames `mat_pval' = $nspec1 $nspec2 $nspec3 $nspec4 
estadd matrix diff=`mat_diff': bdm_tioli
estadd matrix se=`mat_se': bdm_tioli
estadd matrix pval=`mat_pval': bdm_tioli
matrix drop `mat_diff' `mat_se' `mat_pval'

* Create table with esttab
local right_fmt=7*"2 "+1*"%15.0fc "+8*"2 "+	2*"%15.0fc "+4*"2 "	// trick to display no decimal digits for monetary amounts
esttab bdm tioli bdm_tioli using "$results_tab/all_balance_price.tex", replace cells((mean(pattern(1 1 0) fmt(`right_fmt')) diff(pvalue(pval) star)) (sd(pattern(1 1 0) par) se(par))) star(* 0.10 ** 0.05 *** 0.01) ///
	label booktabs nonumbers nomtitles collabels("Mean" "Diff." "Mean" "Diff.") ///
	mgroups("BDM sample" "TIOLI sample" "BDM vs TIOLI", pattern(1 1 1) prefix(\multicolumn{@span}{c}{) suffix(}) span erepeat(\cmidrule(lr){@span}) begin("&(1)&(2)&(3)&(4)&(5)\\")) ///
	substitute("%" "\%" "_" "\_" "$" "\\$")
	
			
*********************************************************************
*** Table 2: EXAMPLES BDM use in developing countries             ***
********************************************************************* 
*IN PDF
	
******************************************************
*** Table 3: Demand at each price                  ***
******************************************************
use "$data_fin/data_anonymized", clear

* Rescale price and relabel original price var
gen price_k=price/1000
lab var price_k "Price (000)"
lab def price 0 "Price = 0" 1000 "Price = 1,000" 3000 "Price = 3,000" 5000 "Price = 5,000" 10000 "Price = 10,000" 15000 "Price = 15,000" 20000 "Price = 20,000", modify

* Adjust full payment indicator
replace paid_full=1 if price==0

* Generate interaction for price and TIOLI indicator
gen tioliprice=price_k*tioli_spl
label var tioliprice "TIOLI * Price"
label var tioli_spl "TIOLI"


*** Demand for BDM sample ***
* Keep only BDM sample
preserve
keep if tioli_spl==0

eststo clear

* Store regression results for willingness to pay and linearity test
eststo wtp1: reg bought_course price_k, rob
eststo wtp2: reg bought_course i.price, rob
test _b[1000.price]=_b[3000.price]/3=_b[5000.price]/5=_b[10000.price]/10=_b[15000.price]/15=_b[20000.price]/20
estadd scalar plin =r(p): wtp2
eststo wtp3: reg bought_course i.price i.form, rob
test _b[1000.price]=_b[3000.price]/3=_b[5000.price]/5=_b[10000.price]/10=_b[15000.price]/15=_b[20000.price]/20
estadd scalar plin=r(p): wtp3
test _b[2.form]=_b[3.form]=_b[4.form]=0
estadd scalar ptreat=r(p): wtp3
summ bought_course if price_k==0, meanonly
estadd scalar ymean=r(mean): wtp*

* Store regression results for paid full price and linearity test
eststo fll1: reg paid_full price_k , rob
eststo fll2: reg paid_full i.price, rob
test _b[1000.price]=_b[3000.price]/3=_b[5000.price]/5=_b[10000.price]/10=_b[15000.price]/15=_b[20000.price]/20
estadd scalar plin =r(p): fll2
eststo fll3: reg paid_full i.price i.form, rob
test _b[1000.price]=_b[3000.price]/3=_b[5000.price]/5=_b[10000.price]/10=_b[15000.price]/15=_b[20000.price]/20
estadd scalar plin=r(p): fll3
test _b[2.form]=_b[3.form]=_b[4.form]=0
estadd scalar ptreat=r(p):  fll3 
summ paid_full if price_k==0, meanonly
estadd scalar ymean=r(mean): fll*

restore

*** Demand for TIOLI sample ***
* Keep only TIOLI sample
preserve
keep if tioli_spl==1

replace form=1 if credit==0
replace form=2 if credit==1

* Store regression results for paid full price and linearity test
eststo tfll1: reg paid_full price_k, rob
eststo tfll2: reg paid_full i.price, rob
eststo tfll3: reg paid_full i.price i.form, rob
summ paid_full if price==0, meanonly
estadd scalar ymean=r(mean): tfll*
restore

* Create table with esttab
esttab _all using "$results_tab/all_demand.tex", replace b(3) se(3) star(* 0.10 ** 0.05 *** 0.01) ///
	stats(N, labels("Observations" ) fmt(%5.0fc  )) ///
	label nocons noomit nobaselevels nonotes nonumbers nomtitles collabels(none) booktabs ///
	mgroups("WTP $\geq$ Price \\ (BDM)" "Paid full price \\ (BDM)" "Paid full price \\ (TIOLI)", pattern(1 0 0 1 0 0 1 0 0 ) ///
		prefix(\multicolumn{@span}{c}{\makecell{) suffix(}}) span begin("&(1)&(2)&(3)&(4)&(5)&(6)&(7)&(8)&(9)\\"))

 
	
*******************************************************
*** Table 4: Correlates of Demand  TIOLI            ***
*******************************************************

use "$data_fin/data_anonymized", clear
keep if tioli_spl==1

*Dummies
replace price=price/1000
lab var price "Price (000)"
gen lprice=log(price)
sum sales_dollars , det
gen lowsales=sales_dollars<r(p50)  
gen highsales=sales_dollars>=r(p50) if sales_dollars!=.
sum risk_taking, det
gen hrisk=(risk_taking>=r(p50)) if risk_taking!=.
gen expect_lowd= (expect_low==1|expect_nochange==1|expect_inc10==1)
sum age, det
gen lowage= (age<r(p50))
gen hage= (age>=r(p50)) if age!=.
sum wealthindex, det
gen pwealth=(wealthindex>=r(p50)) if wealthindex!=.
clonevar m_finance20k_notpos=finance20k_notpos
replace m_finance20k_notpos=0 if missing(finance20k_notpos)
gen d_finance20k_notpos=missing(finance20k_notpos)


cap drop m_* d_*
* Compute variables to replace missing values (keep them in separate globals -e.g. $m_spec1, $d_spec1, etc.)
foreach var in  pwealth highsales edu_comp expect_inc1020 expect_inc20   practice_share pi_index  risk_taking  hage  {
		clonevar m_`var'=`var'
		replace m_`var'=0 if missing(`var')
		gen d_`var'=missing(`var')
	}


label var paid_full "Paid"
label var lprice "Log(price)"
label var m_pwealth "Wealth: above median"
label var m_highsales "Sales: above median"
labe var m_hage "Age: above median"


* Regress WTP on covariates
eststo clear
*  Column 1: gender
eststo: reg paid_full lprice  female   , rob
*  Column 2: wealth
eststo: reg paid_full  lprice  m_pwealth d_pwealth   , rob
* Column 3: sales
eststo: reg paid_full lprice   m_highsales d_highsales  , rob
* Column 4: edu
eststo: reg paid_full lprice   m_edu_comp d_edu_comp , rob
* Column 5: all 4
eststo: reg paid_full  lprice  female m_pwealth m_highsales  m_edu_comp   d_pwealth d_highsales d_edu_comp, rob
* Column 6: all
eststo: reg paid_full lprice  female m_pwealth m_highsales m_edu_comp m_hage m_expect_inc1020 m_expect_inc20   m_practice_share m_pi_index  m_risk_taking d_pwealth d_highsales d_edu_comp  d_hage d_expect_inc1020  d_practice_share d_pi_index  d_risk_taking  , rob
* Column 7: Attributes chosen by lasso
eststo: reg paid_full lprice  female m_pwealth  d_pwealth, rob
estadd ysumm, mean: _all

	  
* Lasso selection
lasso linear paid_full lprice  female m_hage m_pwealth m_highsales m_edu_comp m_expect_inc1020 m_expect_inc20   m_practice_share m_pi_index  m_risk_taking d_pwealth d_highsales d_edu_comp d_hage d_expect_inc1020  d_practice_share d_pi_index  d_risk_taking,  selection(cv, folds(5)) rseed(3214)
lassocoef
*Selected: lprice female m_wealth


* Create table with esttab
esttab _all using "$results_tab/tabf2na", replace b(2) se(2) r2 star(* 0.10 ** 0.05 *** 0.01) drop(d_* _cons) ///
	stats(N r2 ymean, labels("Observations" "R-squared" "Mean dep. var.") fmt(%5.0fc 3 2)) ///
	nobaselevels label booktabs nomtitles collabels(none) nonotes ///
	posthead("&\multicolumn{7}{c}{Paid Full}\\" "\midrule") ///
	substitute("%" "\%" "_" "\_" "$" "\\$")
	
	
*******************************************************
*** Table 5: Correlates of Demand  BDM             ***
*******************************************************

use "$data_fin/data_anonymized", clear
keep if tioli_spl==0

*Dummies
replace price=price/1000
replace wtp_real=wtp_real/1000
lab var price "Price (000)"
gen lprice=log(price)
sum sales_dollars , det
gen lowsales=sales_dollars<r(p50)  
gen highsales=sales_dollars>=r(p50) if sales_dollars!=.
sum risk_taking, det
gen hrisk=(risk_taking>=r(p50)) if risk_taking!=.
gen expect_lowd= (expect_low==1|expect_nochange==1|expect_inc10==1)
sum age, det
gen lowage= (age<r(p50))
gen hage= (age>=r(p50)) if age!=.
sum wealthindex, det
gen pwealth=(wealthindex>=r(p50)) if wealthindex!=.
clonevar m_finance20k_notpos=finance20k_notpos
replace m_finance20k_notpos=0 if missing(finance20k_notpos)
gen d_finance20k_notpos=missing(finance20k_notpos)

*Reneging: defined as not paying full not restricting to buying course, sample is 383 obs.
gen reneged=(paid_full==0) 


cap drop m_* d_*
* Compute variables to replace missing values (keep them in separate globals -e.g. $m_spec1, $d_spec1, etc.)
foreach var in  pwealth highsales edu_comp expect_inc1020 expect_inc20   practice_share pi_index  risk_taking  hage  {
		clonevar m_`var'=`var'
		replace m_`var'=0 if missing(`var')
		gen d_`var'=missing(`var')
	}


label var paid_full "Paid"
label var lprice "Log(price)"
label var m_pwealth "Wealth: above median"
label var m_highsales "Sales: above median"
labe var m_hage "Age: above median"


* Regress WTP on covariates
eststo clear
*  Column 1: gender
eststo: reg paid_full lprice  female   , rob
*  Column 2: wealth
eststo: reg paid_full  lprice  m_pwealth d_pwealth   , rob
* Column 3: sales
eststo: reg paid_full lprice   m_highsales d_highsales  , rob
* Column 4: edu
eststo: reg paid_full lprice   m_edu_comp d_edu_comp , rob
* Column 5: all 4
eststo: reg paid_full  lprice  female m_pwealth m_highsales  m_edu_comp   d_pwealth d_highsales d_edu_comp, rob
* Column 6: all
eststo: reg paid_full lprice  female m_pwealth m_highsales m_edu_comp m_hage m_expect_inc1020 m_expect_inc20   m_practice_share m_pi_index  m_risk_taking d_pwealth d_highsales d_edu_comp  d_hage d_expect_inc1020  d_practice_share d_pi_index  d_risk_taking  , rob
* Column 7: Attributes chosen by lasso
eststo: reg paid_full lprice   m_highsales m_edu_comp    d_pwealth d_highsales  d_edu_comp, rob

*Column 8: WTP
eststo: reg wtp_real   female m_pwealth m_highsales m_edu_comp m_hage m_expect_inc1020 m_expect_inc20   m_practice_share m_pi_index  m_risk_taking d_pwealth d_highsales d_edu_comp  d_hage d_expect_inc1020  d_practice_share d_pi_index  d_risk_taking  , rob
*Column 9: Reneged
eststo: reg reneged lprice  female m_pwealth m_highsales m_edu_comp m_hage m_expect_inc1020 m_expect_inc20   m_practice_share m_pi_index  m_risk_taking d_pwealth d_highsales d_edu_comp  d_hage d_expect_inc1020  d_practice_share d_pi_index d_risk_taking  , rob
estadd ysumm, mean: _all

	 
* Lasso selection for Column 7
lasso linear paid_full lprice  female m_hage m_pwealth m_highsales m_edu_comp m_expect_inc1020 m_expect_inc20   m_practice_share m_pi_index  m_risk_taking d_pwealth d_highsales d_edu_comp d_hage d_expect_inc1020  d_practice_share d_pi_index  d_risk_taking,  selection(cv, folds(5)) rseed(3214)
lassocoef
*Selected: lprice m_highsales m_edu_comp 


* Create table with esttab
esttab _all using "$results_tab/tabf2nb", replace b(2) se(2) r2 star(* 0.10 ** 0.05 *** 0.01) drop(d_* _cons) ///
	stats(N r2 ymean, labels("Observations" "R-squared" "Mean dep. var.") fmt(%5.0fc 3 2)) ///
	nobaselevels label booktabs nomtitles collabels(none) nonotes ///
	posthead("&\multicolumn{7}{c}{Paid Full}" "&\multicolumn{1}{c}{WTP}"  "&\multicolumn{1}{c}{Renege} \\"  "\midrule" )  ///
	substitute("%" "\%" "_" "\_" "$" "\\$")
	


	
*****************************************************************
*** Table 6: Effect of price on attendance                    ***
*****************************************************************

use "$data_fin/data_anonymized", clear

* Rescale price 
gen price_k=price/1000
lab var price_k "Price (000)"

* Relabel original price var
lab def price 0 "Price = 0" 1000 "Price = 1,000" 3000 "Price = 3,000" 5000 "Price = 5,000" 10000 "Price = 10,000" 15000 "Price = 15,000" 20000 "Price = 20,000", modify

eststo clear

keep if price_k!=0 & paid_full==1


*2 outcomes: attending 1 class and classes attended (if paid), for BDM and TIOLI separately
eststo att1: reg atleast1class price_k if tioli_spl==0, rob
eststo att2: reg atleast1class i.price if tioli_spl==0, rob
summ atleast1class if e(sample) & price==1000, meanonly
estadd scalar ymeanb=r(mean): att2
eststo att1b: reg atleast1class price_k if tioli_spl==1, rob
eststo att2b: reg atleast1class i.price if tioli_spl==1, rob
summ atleast1class if e(sample) & price==5000, meanonly
estadd scalar ymeant=r(mean): att2b

eststo att3: reg attendance_total price_k if tioli_spl==0, rob
eststo att4: reg attendance_total i.price if tioli_spl==0, rob
summ attendance_total if e(sample) & price==1000, meanonly
estadd scalar ymeanb=r(mean): att4
eststo att3b: reg attendance_total price_k if tioli_spl==1, rob
eststo att4b: reg attendance_total i.price if tioli_spl==1, rob
summ attendance_total if e(sample) & price==5000, meanonly
estadd scalar ymeant=r(mean): att4b


* Create table with esttab
esttab att* using "$results_tab/attendancebis", replace b(3) se(3) star(* 0.10 ** 0.05 *** 0.01) drop(_cons) ///
	stats(N ymeanb ymeant, labels("Observations"  "Mean Y for P=1,000"  "Mean Y for P=5,000") fmt(%5.0fc 3 3)) ///
	nobaselevels label booktabs nomtitles collabels(none) nonotes ///
	posthead("&\multicolumn{4}{c}{At least 1 class}" "&\multicolumn{4}{c}{N. classes}\\"  "&\multicolumn{2}{c}{(BDM)}" "&\multicolumn{2}{c}{(TIOLI)}" "&\multicolumn{2}{c}{(BDM)}" "&\multicolumn{2}{c}{(TIOLI)} \\"  "\midrule" )  ///
	substitute("%" "\%" "_" "\_" "$" "\\$") 

 			
*********************************************************************
*** Table 7: Cost-Effectiveness                                   ***
********************************************************************* 
*IN PDF
 


*****************************************
*** Table A1:  Statistics on Payments ***
*****************************************

use "$data_fin/data_anonymized", clear

* Generate copies of variables, with missing values (for conditional subsamples)
clonevar attend_tot_ifpaid=attendance_total if paid_full==1
lab var attend_tot_ifpaid "N. of classes (if paid full)"

* Summary stats for BDM sample
eststo clear
foreach p in 0 1000 3000 5000 10000 15000 20000 {
	estpost tabstat  paid_something paid_full  atleast1class  atleast5class ///
		allclasses  attendance_total  attend_tot_ifpaid if tioli_spl==0 & price==`p', column(stat)
	eststo bdm_`p', title("`: dis %9.0fc `p''")		// trick to save the (formatted) price as title and then use it with esttab
}
estpost tabstat  paid_something paid_full  atleast1class  atleast5class ///
	allclasses  attendance_total  attend_tot_ifpaid if tioli_spl==0, stat(count) column(stat)
eststo bdm_obs, title("Obs.")		// trick to set the title to use with esttab

* Summary stats for TIOLI sample
foreach p in 0 5000 10000 15000 {
	estpost tabstat paid_something paid_full  atleast1class  atleast5class ///
		allclasses attendance_total  attend_tot_ifpaid if tioli_spl==1 & price==`p', column(stat)
	eststo tio_`p', title("`: dis %9.0fc `p''")		// trick to save the (formatted) price as title and then use it with esttab
}
estpost tabstat paid_something paid_full  atleast1class  atleast5class ///
	allclasses attendance_total attend_tot_ifpaid if tioli_spl==1, stat(count) column(stat)
eststo tio_obs, title("Obs.")		// trick to set the title to use with esttab

* Create table with esttab
esttab bdm* using "$results_tab/all_summary_byprice_transpose", replace cells((mean(pattern(1 1 1 1 1 1 1 0) fmt(2)) count(pattern(0 0 0 0 0 0 0 1) fmt(%10.0fc)))) ///
	label booktabs nonumbers mtitles collabels(none) ///
	mgroups("Price offered", pattern(1 0 0 0 0 0 0 1) prefix(\multicolumn{@span}{c}{) suffix(}) span begin("&(1)&(2)&(3)&(4)&(5)&(6)&(7)&(8)\\") end(\cmidrule(lr){2-8})) ///
	posthead("\midrule" "\multicolumn{9}{l}{\textbf{Panel A. BDM sample}}\\") postfoot("")
esttab tio* using "$results_tab/all_summary_byprice_transpose", append cells((mean(pattern(1 1 1 1 0) fmt(2)) count(pattern(0 0 0 0 1) fmt(%10.0fc)))) ///
	extracols(2 2 5) label booktabs nonumbers nomtitles collabels(none) ///
	prehead("\midrule") posthead("\multicolumn{9}{l}{\textbf{Panel B. TIOLI sample}}\\")

	
*********************************************************************
*** Table A2: IV Effect on main outcomes, BDM sample              ***
********************************************************************* 

use "$data_fin/data_anonymized", clear
cap ssc install ivreg2
cap ssc install ranktest

keep if tioli_spl==0

gen price_k=price/1000
lab var price_k "Price (in J\\\$000s)"
lab def price 0 "Price = 0" 1000 "Price = 1,000" 3000 "Price = 3,000" 5000 "Price = 5,000" 10000 "Price = 10,000" 15000 "Price = 15,000" 20000 "Price = 20,000", modify

* Compute variables to replace missing values (covariates for IPW calculation)
global vars_ipw
foreach var in $spec1 $spec2 $spec3 $spec4 sales_lm profits_lm {
	cap confirm new variable m_`var'
	if _rc continue
	clonevar m_`var'=`var'
	replace m_`var'=0 if missing(`var')
	gen d_`var'=missing(`var')
	global vars_ipw "$vars_ipw m_`var' d_`var'"
}

* Compute the inverse probability weights
gen f_natt=1-f_attriter
probit f_natt $vars_ipw 
predict pscore

gen ipw=1/pscore
lab var ipw "Inverse probability of attrition weight"
drop f_natt pscore*

* Define WTP groups in one categorical variable
egen wtp_group=cut(wtp_real), at(0, 3001,5001,10001,15001,20001,1000000) icodes label
lab def wtp_group ///
	0 "0 \$\leq\$ WTP \$\leq\$ 3,000" ///
	1 "3,000 \$<\$ WTP \$\leq\$ 5,000" ///
	2 "5,000 \$<\$ WTP \$\leq\$ 10,000" ///
	3 "10,000 \$<\$ WTP \$\leq\$ 15,000" ///
	4 "15,000 \$<\$ WTP \$\leq\$ 20,000" ///
	5 "WTP \$>$ 20,000", modify
lab val wtp_group wtp_group
* Generate dummies for WTP groups
for num 1/5: gen wtpgX=wtp_group==X


* Run regressions for main outcomes with IPW
eststo clear
foreach var in practice_share pi_index knowledge has_employees sales_pos_lm profits_pos_lm {
	eststo ipw_`var': ivreg2 f_`var' (atleast1class=price_k) wtpg1-wtpg5 m_`var' d_`var' [weight=ipw], rob
}
foreach var in win_sales_lm win_profits_lm {
	eststo ipw_`var': ivreg2 f_`var' (atleast1class=price_k) wtpg1-wtpg5 m_`var' d_`var' [weight=ipw], rob
}


*Interactions
sum wtp_real if f_knowledge!=. & ipw!=.
gen wtp_demean=(wtp_real-r(mean))/1000
* Generate interactions
gen atleast1class_wtp=atleast1class*wtp_demean
gen price_wtp=price*wtp_demean
label var atleast1class_wtp "At least 1 class * WTP demeaned"
label var price_wtp "Price * WTP demeaned"

* Run regressions for interactions with IPW
foreach var in practice_share pi_index knowledge has_employees sales_pos_lm profits_pos_lm {
	eststo nipw_`var': ivreg2 f_`var' (atleast1class atleast1class_wtp =price price_wtp) wtp_demean wtpg1-wtpg5 m_`var' d_`var' [weight=ipw], rob
	summ f_`var' if e(sample) & atleast1class==0
	estadd scalar ymean1=r(mean)
}
foreach var in win_sales_lm win_profits_lm {
	eststo nipw_`var': ivreg2 f_`var' (atleast1class atleast1class_wtp =price price_wtp) wtp_demean wtpg1-wtpg5 m_`var' d_`var' [weight=ipw], rob
	summ f_`var' if e(sample) & atleast1class==0
	estadd scalar ymean2=r(mean)
}


* Create table with esttab
esttab ipw* using "$results_tab/bdm_iv.tex", replace b(2) se(2) star(* 0.10 ** 0.05 *** 0.01) ///
	cells((b(pattern(1 1 1 1 1 1 0 0) fmt(2) star) b(pattern(0 0 0 0 0 0 1 1) fmt(%10.0fc) star)) ///
		(se(pattern(1 1 1 1 1 1 0 0) fmt(2) par) se(pattern(0 0 0 0 0 0 1 1) fmt(%10.0fc) par))) ///
	keep(atleast1class) nobaselevels label booktabs collabels(none) nonotes ///
	stats(r2, labels("R-squared") layout(@) fmt(3)) ///
	mlabels("Business\\ practices" "Personal \\ initiative" "N. correct \\ answers" "Has \\ empl." "Pos. \\ sales" "Pos. \\ profits" "Sales" "Profits", prefix(\makecell{) suffix(})) ///
	posthead("\midrule" "\multicolumn{9}{l}{\textbf{Panel A. Main Effects}}\\" "\addlinespace") postfoot("") ///
	substitute("_" "\_" "$" "\$" "@" "")
esttab nipw* using "$results_tab/bdm_iv.tex", append b(2) se(2) star(* 0.10 ** 0.05 *** 0.01) ///
	cells((b(pattern(1 1 1 1 1 1 0 0) fmt(2) star) b(pattern(0 0 0 0 0 0 1 1) fmt(%10.0fc) star)) ///
		(se(pattern(1 1 1 1 1 1 0 0) fmt(2) par) se(pattern(0 0 0 0 0 0 1 1) fmt(%10.0fc) par))) ///
	stats(N r2 ymean1 ymean2, labels("Observations" "R-squared" "Mean Y if did not attend any class") layout(@ @ "@@") fmt(%5.0fc 3 2 %10.0fc)) ///
	keep(atleast1class atleast1class_wtp ) nobaselevels label booktabs nonumbers nomtitles collabels(none) nonotes ///
	prehead("") posthead("\midrule" "\multicolumn{9}{l}{\textbf{Panel B. Heterogeneous effects by WTP}}\\" "\addlinespace") ///
	substitute("_" "\_" "$" "\$" "@" "")



*MULTIPLE HYPOTHESES CORRECTION*
** Wyoung controlling for baseline values
cap scc install wyoung
wyoung, cmd("ivreg f_practice_share  (atleast1class=price_k) wtpg1-wtpg5 m_practice_share d_practice_share [weight=ipw], rob" ///
"ivreg f_pi_index (atleast1class=price_k) wtpg1-wtpg5 m_pi_index d_pi_index [weight=ipw], rob" ///
"ivreg f_knowledge (atleast1class=price_k) wtpg1-wtpg5 m_knowledge  d_knowledge  [weight=ipw], rob" ///
"ivreg f_has_employees (atleast1class=price_k) wtpg1-wtpg5 m_has_employees d_has_employees [weight=ipw], rob" ///
"ivreg  f_sales_pos_lm (atleast1class=price_k) wtpg1-wtpg5 m_sales_pos_lm d_sales_pos_lm [weight=ipw], rob" ///
"ivreg    f_profits_pos_lm (atleast1class=price_k) wtpg1-wtpg5 m_profits_pos_lm d_profits_pos_lm [weight=ipw], rob" ///
"ivreg f_win_sales_lm  (atleast1class=price_k) wtpg1-wtpg5 m_win_sales_lm  d_win_sales_lm  [weight=ipw], rob" ///
"ivreg f_win_profits_lm (atleast1class=price_k) wtpg1-wtpg5 m_win_profits_lm d_win_profits_lm [weight=ipw], rob") ///
familyp(atleast1class) bootstraps(1000) seed(20)

** Wyoung controlling for baseline values
wyoung, cmd("ivreg f_practice_share  (atleast1class atleast1class_wtp =price price_wtp) wtp_demean wtpg1-wtpg5 m_practice_share d_practice_share [weight=ipw], rob" ///
"ivreg f_pi_index (atleast1class atleast1class_wtp =price price_wtp) wtp_demean wtpg1-wtpg5 m_pi_index d_pi_index [weight=ipw], rob" ///
"ivreg f_knowledge (atleast1class atleast1class_wtp =price price_wtp) wtp_demean wtpg1-wtpg5 m_knowledge  d_knowledge  [weight=ipw], rob" ///
"ivreg f_has_employees (atleast1class atleast1class_wtp =price price_wtp) wtp_demean wtpg1-wtpg5 m_has_employees d_has_employees [weight=ipw], rob" ///
"ivreg  f_sales_pos_lm (atleast1class atleast1class_wtp =price price_wtp) wtp_demean wtpg1-wtpg5 m_sales_pos_lm d_sales_pos_lm [weight=ipw], rob" ///
"ivreg    f_profits_pos_lm (atleast1class atleast1class_wtp =price price_wtp) wtp_demean wtpg1-wtpg5 m_profits_pos_lm d_profits_pos_lm [weight=ipw], rob" ///
"ivreg f_win_sales_lm  (atleast1class atleast1class_wtp =price price_wtp) wtp_demean wtpg1-wtpg5 m_win_sales_lm  d_win_sales_lm  [weight=ipw], rob" ///
"ivreg f_win_profits_lm (atleast1class atleast1class_wtp =price price_wtp) wtp_demean wtpg1-wtpg5 m_win_profits_lm d_win_profits_lm [weight=ipw], rob") ///
familyp(atleast1class_wtp) bootstraps(1000) seed(20)	
	
	




